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A curved bar structural analysis is developed for an rc-degree elliptical elastic ring with sinusoid- 
ally varying cross-section dimensions, where the n-degree ellipse is denned by 
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1, /i 2*2.0. 



A sequential grid constrained optimization method is used to search for a minimum weight design in 
six dimensional shape parameter space. Numerical results indicate that rings of this shape can be 
designed to have considerably less weight and greater flexibility than comparable circular rings with 
uniform cross-section dimensions. 
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programming, optimization, proving ring, structural design, theory. 



1. Introduction 

Elastic ring force transducers (such as proving rings, load rings, ring dynamometers, etc.), 
for a wide range of load carrying capacities have most frequently been made circular in shape with 
uniform rectangular cross section. This shape has imposed limitations on the usefulness of this 
type device, particularly for higher capacities (i.e., greater than 100,000 lbf), due to overall size, 
gross weight, limited deflection range, and size effect during heat treatment. For the same force 
resisting capacity, rings of more complex shapes can be designed to have greater deflection ranges, 
less weight, and smaller maximum dimensions. This paper describes methods of analysis of a 
broad geometric class of complex shaped rings which result in improved deflection and weight 
characteristics. 

Structural analysis and design optimization methods are developed for an Ai-degree elliptical 
elastic ring with sinusoidally varying cross-section dimensions. The ring is analyzed as a thick 
curved bar, assuming small deflections, plane strain, Winkler-Bach bending stress distribution, 
uniform normal stress distribution, and parabolic shear stress distribution [l]. 1 A sequential grid 
constrained optimization method is used to search for the minimum weight ring design that will 
satisfy the specific force capacity, deflection, and dimensional requirements of a force transducer. 



1 Figures in brackets indicate the literature references on page 150. 
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List of Symbols 

x, y rectangular coordinates of ellipse, 

r, 6 polar coordinates of ellipse, 

n degree of ellipse, 

a, b semimajor and semiminor axes of ellipse, 

h ring cross-section width, 

t ring cross-section thickness normal to tangent of ellipse, 

c average of maximum and minimum ring width, 

d average of maximum and minimum ring thickness, 

e dimensionless parameter of ring width variation, 

/ dimensionless parameter of ring thickness variation, 

5 arc length variable measured along ellipse, 

i// counterclockwise angle from radius vector to tangent line, 

L total load applied to ring, 

P = L/2, load applied to one quadrant of ring, 

M moment, 

TV normal force, 

V shearing force, 

a normal stress, 

r shearing stress, 

e strain, 

k curvature of ellipse, 

/3 relative rotation of differential segment ends due to TV, 

G curved bar rotation, 

8 curved bar deflection, 

A total ring deflection, 

U strain energy of curved bar, 

u strain energy per unit volume, 

v Poisson's ratio, 

E modulus of elasticity, 

g integral defined by eq (19), 

S design maximum stress, 

A scale factor, 

W ring weight, 

y material weight per unit volume, 

gi constraint function, 

bi constraint bound, 

Xj variable of the optimization problem, 

z objective function. 

2. Structural Analysis 

An elastic ring force transducer is loaded by two forces acting in opposite directions along 
a diameter, as shown in figure 1. Corresponding ring deflections are measured and related to 
applied forces by a calibration factor. 2 For a ring of the geometric class represented in figure 
1, the locus of centroids of ring cross sections is an n-degree ellipse defined by the equation 

jtf+GT-i. »>2.o, (1) 



2 Wilson, Tate, and Borkowski [2] have described the uniform cross section, circular force transducer, its calibration, and performance under various conditions 
of use. 
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FIGURE 1. 1,000,000 Ibf capacity n-degree elliptical ring with sinusoidally 
varying thickness. 



or in polar coordinate form 



r=ab(b" cos" + a"sin n 6)~ lln . 



(2) 



Several n-degree ellipses are plotted in figure 2. Ring rectangular cross-section width h and thick- 
ness t are given by 

h = c(l— ecos 20), (3) 



t = d(l-fcos 20), 



(4) 



in which c and d are positive parameters and e and /are parameters of absolute value less than 
unity. Using these equations, the shape of a ring within this geometric class can be specified by 
the seven parameters a, b, c, d, e, f, and n. If a = b, n = 2, and e =f= 0, the ring shape is circular 
with uniform rectangular cross section. 

One quadrant of the loaded ring is shown in figure 3. Due to symmetry there is no rotation at 
= or at = 77/2, no shear force Fat 6 = 0, and no normal force /V at 6 = 7t/2. Expressions for the 




FIGURE 2. n-degree elliptical curves. 




FIGURE 3. Loads acting on one quadrant of a ring. 
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resultant moment M and forces N and V acting on a typical normal cross section are derived by 
statics using formulas from calculus [3] as follows: 



dO _ sin \\f 
ds r 



1 



-'a)" 



(5) 



dr _ _dr dO 

ds dO ds ' 



M = M„i2 - Pr cos 9 = PM' , 
N= P sin (6 + xl>) = P (sinOj+ r cos 8 j) = PN' , 



(6) 
(7) 
(8) 



rsinOj--cosej)=PV. 



Differentiation of eq (2) gives 



dr_ (b n cos*""" 6 sin 6 -a" sin'"-'> 6 cos g) 
dd~ r (6" cos" 0+a"sin"0) 



(9) 



(10) 



The assumption that ring cross sections remain plane leads to the Winkler-Bach bending stress 
formula, which has been written by Seely and Smith [1] in the following form, valid for a straight 
or curved bar: 



_M 



: + - 



y 



(l + ky) 



/.«! 



+ ky) 



dy 



(id 



The variable y is the distance from the cross-section centroid to a point in the cross section, positive 
in the outward direction. The curvature k, the inverse of the radius of curvature, is given by [3] 



'■+2 ra 



k=- 



dr\ z d 2 r 



del do 1 



[-(£)'] 



3/2 



(12) 



Differentiation of eq (10) gives the required second derivative 



—;=(n+l)ab[a n sin*"- 1 ) cos 0-6" cos*"- 1 ) 6 sin 6] 2 [a" sin" 6 + b n cos" 0]-<2 + i/*) 

+ « ( " + 1) 6[sin" 6-(n-\) sin ( "- 2) cos 2 0] [a" sin" + 6" cos" 0]-< 1 + 1 /") 
+ a6 ( " + 1) [cos" O-(n-l) cos ( "- 2) sin 2 0] [a" sin" (9+6" cos" 0]- (1 + 1 /"). 



(13) 



It is assumed that stress due to the normal force N is uniformly distributed over the cross 
section, that is 

N 



<Jn = 



ht 

142 



(14) 



It is assumed that shear stress is distributed parabolically over the rectangular cross section 
according to the formula 



V /3 



T V ~- 



ht \2 



~6 



r 



(15) 



Castigliano's strain energy theorem [1] is used to determine the ring bending moment and 
deflection at d—ir\2. The theorem stated symbolically as applied to this problem is 



077-/2 : 

0\r/2 : 



dU 

dU 

dP' 



(16) 



(17) 



in which U = total strain energy of the loaded bar, 

Oir/2 = bar rotation at 7r/2, 
and 877-/2 — bar deflection at tt/2 in direction of P. 

The strain energy per unit volume due to bending stress is 



1 

Um — ~Z 0~a/€a/ : 



(l—V 2 )(TJ f 

2E 



(18) 



in which €m = unit plane strain, 

^=Poisson's ratio, 
and E= modulus of elasticity. 

Substituting eq (11) into (18) and integrating over a differential segment gives the corresponding 
strain energy of the segment 



dU M = 



2Eh 



k 



■+■ 



y 



tt + ky) 






ky) 



dy 



'- ) (\-v 2 )M 2 . 

<lv \' h ~ ~2Eh — gds - 



The strain energy per unit volume due to the normal force N is 



(19) 



_1 _1 (N\ ( N \_ N 2 

u N -- <r N € N - ^ 1^1 1^1 - ^^, 



(20) 



and the strain energy of a differential segment is 



dU " = 2Eht ds - 



(21) 



The strain energy per unit volume due to the shear force V is 

(1 + ^t 2 (l + i>)F/3 yV 
0-672 . 



Uv 



E EhH 2 \2 " t 2 ) 

and the strain energy of a differential segment is 



6 /l + v\ V 2 



5\ E } ht 
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(22) 



(23) 



The normal force N acting on a differential segment causes rotation of one end cross section 
relative to the other, due to the difference in arc length from inside radius to outside radius. The 
mean change in segment arc length is 



Eht \k 



-G 



(i- 



(24) 



in which /3 = rotation due to N. 

The bending moment M acting during this rotation contributes strain energy to the segment 
equal to 



dU MN = - 



kMNds 
ThT' 



(25) 



Adding eqs (19), (21), (23), and (25) gives the total strain energy of a differential segment 



1 rA r2 6 PI M 2 g MNk 



ds. 



(26) 



Applying eq (16) and the boundary condition of no rotation at 6 = ttJ2 gives 



du = i r 12 

dMrrlt E Jo 



TV dN 



12 



ht dM n /2 5 



+ -E-(l + y) 



V 3V 



ht dMirl 



•+(1-^)-^ 



gM dM kM dN kN dM 



•+ 



+- 



h dM*/* ht dM n!2 ht dM % 



ds 



\: 



1 fl- \(\-v-)gM kN' 
h + JT 



ds = 0. 



(27) 



Substitution of eqs (7) and (8) into (27) gives 

(1 — v-)gr cos 6 



L 



M*, t = P 



Applying eq (17) to (26) gives 



f k sin ( k, 

ds -J -J^ dr -\ ~ 



r cos 6 



ht 



de 



/ 



(i - p 2 )g 



ds 



(28) 



du i r n 



_du_ j_ r 
■' 2 ~Hp~e)„ 



N dN l2(l + v)V dV 



ht flP. + " 



gM dM kM dN kN dM 



E J,, ht 



5ht dP + {1 vl) h 3P + ht c)P + ht dP 



12 



ds 



W' 2 + y (l + v)V' 2 - (l -u 2 )gtr cos dM' + kM'N'-kr cos ON' 



ds. 



(29) 



The tangential stress on the inner surface is provisionally assumed to be the limiting stress 
that determines the force resisting capacity of the ring. The possibility of capacity being limited 
by some other component of stress or combined stress should be checked in any final design 
analysis. Substituting — 1\1 for y outside the integral in eq (11) and adding to eq (14) gives the 
tangential stress on the inner surface of the ring 



a i 



M'k 
hT 



M't 



h(2-kt) 



-I 



y 2 

(1 + ky) 



dy 



N' 
+ Yt 



(30) 
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For a circular ring of uniform cross section the maximum value of cry occurs at = 7r/2. But for the 
more general geometry considered here the maximum may occur at = or at some point in 
between. 

The structural analysis thus far has been developed for one quadrant of the ring. For the total 
ring, figure 1, the applied force is 

L=2P, (31) 

and the deflection at the point of load application is 

A = 28^/2. (32) 

3. Design Analysis 

The design analysis attempts to find the dimensions of the minimum weight ring that will 
satisfy the arbitrary force capacity, deflection, and dimensional requirements of a particular force 
transducer. The deflection range of the ring must be compatible with the deflection sensing device, 
and the inside dimensions of the ring must be sufficient to accommodate the sensor. Outside di- 
mensions of the ring must not interfere with its intended use. The thickness of the ring must permit 
the heat treatment required for the ring material. All trial solutions are scaled up or down to 
coincide with the force capacity requirement. 

3.1. Scale Factors 

The procedure adopted here for meeting the force capacity requirement is to scale all length 
dimensions in the r-0 plane by a factor A that will make the maximum stress (the maximum value 
of &t for the entire ring) equal a prescribed maximum design stress S. The cross-section width h 
is not scaled. This gives the specialized equation for the maximum tangential stress on the inner 
surface of the scaled ring 



* 2 



G) 


'M'k 
ht 


UF 




2X ' 





M't at 

ht 



w-^La+P 



ky) 



(33) 



in which S = design maximum stress, 

L\ = design capacity, 
and F= maximum value of bracketed function. 

This gives the scale factor 

Scaled ring dimensions, A subscripted, are 

r\=kr, (35) 

and 

tx=kt. (36) 
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Scaled forces and moment are 



Scaled ring deflection is 



Scaled ring weight is 



U = kL, (37) 

Nx = \N, (38) 

Vx = \V, (39) 

Mx = \ 2 M. (40) 

Ax = AA. (4i) 

72 



W k = ^y\ 2 htds=kW (42) 

in which y= material weight per unit volume. 

3.2. Thin Circular Ring 

By making the appropriate simplifying assumptions, the thick ring equations developed above 
can be reduced to relatively simple thin circular ring equations. The thin circular ring equations, 
although they may be in error by several percent, can be used to determine a reasonable range of 
ring proportions for beginning an optimum search procedure. The equations also give a useful 
indication of the comparative efficiencies of different ring materials. 

A thin ring is here defined as one in which the ratio of thickness to radius of curvature is so 
small that, with negligible error: (1) the stress distribution due to bending can be assumed linear 
over a cross section, and (2) deflection due to shear V and normal force TV can be neglected. By 
assuming zero effective curvature and no deflection due to V or N, the following thin ring equations 
for moment, maximum bending stress, and deflection at the point of load application can be 
obtained from eqs (28) through (32): 

M = -, (43) 

and 



These equations, solved simultaneously, give 



24\ (l~-v%r* 

Eht* ' (45) 



3/ 72EL 2 A ,, 

V(7r 4 -877 2 )(l-^)c^ 3 7 ^ 2, {qf ° } 



and 

7ro-r/^ 2 



r = 



6L 
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(47) 



The weight of a thin circular ring is 



24£yLA 



r = 2 ^ ^-8 )( WK m 

Note that the weight, for a particular LA, is a function of the material properties E, v, cry, and y. 
Equation (48) suggests the importance of maximum permissible design stress in a comparison of 
different ring materials where weight is a significant factor. 

3.3. Optimization Method 

The optimization problem considered here is of the general nonlinear programming class. 
The goal of this type problem is to determine the set of values of the variables X\, . . ., x p which 
satisfies the constraint conditions 

gi(x u • • ., x„){ ^ , = , ^ }bt, i= 1, . . ., m, (49) 

and minimizes the objective function 

z=f(x u . . .,%). (50) 

In the present case the p variables are the six independent ring shape parameters a, c, d, e, 
/, and n. The structural analysis is first done for a ring of unit value of mean radius, that is 

a + b 

-2-=l- (51) 

Thus the value of the parameter a determines the value of b. The thickness parameter d is, for the 
unit size ring, the ratio of the mean of the extreme values of thickness to mean radius. After the 
structural equations are solved for the unit mean radius ring (using Simpson numerical integration), 
the results are scaled to full size using the scale factor X. The width parameter c is not scaled 
during the entire analysis, and the shape parameters e, /, and n are nondimensional and are not 
scaled. 

The m constraint functions g\ and the corresponding constraint bounds b\ are as follows: 

(1) Ring deflection 

gl = Ax^b u (52) 

(2) Maximum thickness of ring cross section 

g 2 = Xrf(l+l/1) ^62, (53) 

(3) Outside width of ring 

g., = k(2a + d-df) ^63, (54) 

(4) Outside height of ring 

# 4 = M4-2a + d+ df) ^ b 4 , (55) 
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(5) Maximum width of ring cross section 

£5=c(l + H)^fc 5 , (56) 

(6) Inside clear height of ring 

g« = k(4-2a-d-df)^b 6 . (57) 

The constraint bounds bi are the limits on ring characteristics imposed by use, fabrication, or 
material requirements. 

The objective function z to be minimized is the ring weight W\, eq (42). 

A sequential grid method was used with reasonable success to search for minimum weight 
solutions. The method requires several computer runs of moderate duration. Between computer 
runs changes are made in the input data based on the experience of previous runs. The "Elastic 
Ring Optimization Program" listed in the appendix was used for this search procedure. 

For the first computer run, one or more values of each of the six variables Xj are chosen as 
input data. If two or more values of a variable are chosen, they are widely spaced, but within a 
reasonable range for that variable. By nested do loops the computer generates all combinations of 
the variables. Each combination is a solution point in the six dimensional shape parameter space. 
The ring weight and the six constraint function values are computed and printed for each solution 
point. The best solution point for a computer run is the one with the least weight that also falls 
within all six constraint bounds. 

Subsequent computer runs are made using as input data variable values distributed about the 
best known solution point. The range of values of any variable should include the value of that 
variable for the best known solution point. Any results of previous computer runs may be useful 
as a guide for selecting input data. The process is repeated until a solution point with satisfactorily 
small weight is found, or until results indicate that a satisfactory approximation of a minimum 
weight solution has been found. 

Computer time required for each solution point is approximately one-half second on a UNIVAC 
1108 computer. Obviously there are practical limits to the number of variable combinations one 
would want to generate in one computer run. 

A second computer program (the "Detailed Ring Analysis Program" listed in the appendix) 
can be used to compute scaled dimensions, forces, moments, and stresses for a particular set of 
input shape parameters. The program prints r\, t\, h, k, N\, V\, M\, err, and the value of the integral 
in eq (11) for each incremental value of 6. These results can be used with eqs (11), (14), and (15) to 
determine whether ring capacity is limited by a stress condition other than the maximum value 

of <Tt> 

The following numerical example demonstrates the optimization method. 

3.4. Numerical Example 

The problem is to determine, for a ring of 1,000,000 lbf capacity and for the material constants 
and constraint bounds listed below: 

(1) The approximate weight of a circular ring of uniform rectangular cross section. 

(2) The dimensions of an approximately minimum weight n-degree elliptical ring of constant 
cross section width and sinusoidally varying thickness. 

(3) The dimensions of an approximately minimum weight rc-degree elliptical ring of sinusoi- 
dally varying cross section width and thickness. 

The material constants are: S= 150,000 lb/in 2 , £ = 30,000,000 lb/in 2 , */=0.3, andy = 0.29 lb/in 3 . 
The constraint bounds are: b\ = 0.25 in, b 2 = 4.5 in, 6 3 = 40.0 in, b 4 = 40.0 in, b*> = 14.0 in, be = 9.0 in. 

a. Solution 1 

Thin circular ring equations were used to compute a rough approximation of ring weight and 
dimensions. Substitution of L, £\ y, v, S for err, and b\ for A into eq (48) gives an approximate ring 

148 



weight of 1365 lb. Substitution of the same values, along with />.-> for //, into eqs (46) and (47) give 
£ = 3.66 in and r= 14.73 in. The ratio of thickness to radius is 0.248. The more accurate thick 
curved bar equations were used for all subsequent computations. 

Computer run No. 1 of the sequential grid search included nine uniform cross-section circular 
rings with widths of 10, 12, and 14 in and thickness to radius ratios of 0.20, 0.2S, and 0.30. For 
these nine cases the ratios of weight to deflection differed by less than two percent, a result that 
might have been anticipated from eq (48). The two cases with characteristics nearest the constraint 
bounds had deflections of 0.234 and 0.277 in and weights of 1327 and 1587 lb respectively. Linear 
interpolation between these two cases gives an estimated weight of 1424 lb for a ring that would 
deflect 0.250 in. 

b. Solution 2 

A sequence of five computer runs were used to seek the dimensions of a minimum weight, 
constant width ring. The input shape parameter data and the resulting weight for the best point of 
each run are given in table 1. The best point parameter values for each run are underlined. The 
weight of 553 lb for the fifth run is about 39 percent of the estimated weight of 1424 lb for a uniform 
cross-section circular ring. The 553 lb ring was accepted as a satisfactory approximation of the 
minimum weight constant width ring. The ring shape is shown in figure 1. Critical ring dimensions 
and a deflection of 0.257 in all fall within the constraint bounds. 

Plots of cry versus 6 for the 553 lb ring and for a 1587 lb uniform circular ring (c=14.0 in, 
<7 = 0.25 in, A = 0.277 in in computer run No. 1) are shown in figure 4. The high stress level over a 
greater portion of the 553 lb ring is an indication of more efficient use of material in this more 
complex ring shape. 

c. Solution 3 

Computer runs 6 and 7 included many cases with variation of cross-section width. But for both 
runs the best points were constant width cases. The seventh run included a variable width case 
with parameter values c= 13.6 in, e = 0.025, and/=0.45, and a resulting weight of 539.18 lb which 

TABLE 1. Sequential grid search data for numerical example 







Independent ring * 


hape parameters* 




Weighl 
for besl 


run No. 


a 


c 


d 


e 


/ 


n 


point 


1 


Inch 
1. 


Inch 
10. 


Inch 
.2 


0. 


0. 


2. 


l.h 
802.00 




1.2 


12. 
14. 


.23 
.3 




2 


3. 
4. 




2 


1.1 


13. 


.25 


0. 


.1 


2.5 


828.42 




1.3 


14. 


.3 
^35 




.2 
.3 


3. 
3.5 




3 


1.15 

1.25 


13.5 

14. 


.28 
.3 


0. 


.2 
.3 


3. 
3.5 


639.84 








.32 




.4 


4. 




4 


1.23 


13.4 


.31 


0. 


.3 


4. 


572.90 




1.25 


13.7 


.33 




.35 


4.5 






1.27 








.4 


5. 




5 


1.22 
1.23 


13.6 

13.8 


.303 
.31 


0. 


.4 

.45 


4.5 


553.06 




1.24 




.315 






_ 




6 


1.21 


12. 


.3 


-.1 


.45 


5.5 


564.55 




1.22 


12.7 


.305 


0. 
.1 


•1 


6. 




7 


1.21 


12.7 

13. 

13.3 

13.6 

14. 


.305 


- .025 

0. 

~025 
.05 
.075 
.1 


.4 
.45 


6. 


539.08 



UNIFORM CIRCULAR RING . 




NONUNIFORM n-DEGREE ELLIPTICAL RING 
(n=5.5) 



10 20 30 40 50 60 70 80 90 

ANGULAR LOCATION, 9 , DEG 

FIGURE 4. Tangential stress on inner surface of 
nonuniform n-degree elliptical ring and uniform 
circular ring. 



1 Besl point values for each computer run arc underlined. 



292-251 0-68— 3 
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is only slightly greater than the best point weight. Additional searching might result in a slight 
weight reduction due to width variation. But the weight savings would probably not be enough 
to justify the additional complexity of shape. 

4. Summary 

A structural analysis method has been developed for an ^-degree elliptical elastic ring with 
sinusoidally varying cross-section dimensions. A sequential grid constrained optimization method 
has been used to search for the minimum weight ring of this geometric class that would be suitable 
as the load supporting element of a force transducer. 

Limited numerical results were obtained for one-million pound force capacity rc-degree ellip- 
tical rings with nonuniform cross-section dimensions. The results indicate that such a ring with 
constant cross-section width and sinusoidally varying thickness can be designed to weigh as little 
as 39 percent of the weight of a comparable uniform cross section, circular ring. Outside dimensions 
of the more complex ring are also significantly less than those of the uniform, circular ring. The 
results for an n-degree elliptical ring with sinusoidally varying cross-section width and thickness 
suggest that the weight savings due to variation of width may not be great enough to justify that 
additional complexity of shape. 

The Elastic Ring Optimization Program (given in appendix) can be used directly for design 
optimization by the sequential grid method described in this paper; or the program could be used 
in combination with some other optimization procedure. In either case, the Detailed Ring Analysis 
Program (see appendix) or another comprehensive stress analysis should be used for a detailed 
examination of stresses for the entire ring. 



The analysis was coded in FORTRAN V for the UNI VAC 1108 computer by R. M. Slesser 
of the National Bureau of Standards Center for Computer Sciences and Technology. 
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6. Appendix 

The two FORTRAN V computer programs used in the elastic ring design analysis are listed 
below. The Elastic Ring Optimization Program is used to search for the optimum set of shape 
parameters a, c, d, e, /, and n. The Detailed Ring Analysis Program is used to compute scaled 
dimensions, forces, moments, and stresses for a particular set of input shape parameters. Table 
Al lists the computer symbols used to represent mathematical symbols appearing in the equations 
of the text. Sample input data is included at the end of each program. 
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FABLE Al. Corresponding symbols in text, computer programs, and 
computer output 



Text 


Programs 


Output 


a 


CI 


A 


c 


C3 


C 


d 


C5 


D 


e 


C7 


E 


f 


C8 


F 


n 


XN 


N 


L k 


XLLAM 


L LAMBDA 


S 


S 


S 


E 


E 


E 


v 


XNHU 


NU 


y 


GAMMA 


GAMMA 


h 


B(U) 


L 


t 


T(IJ) 


T 


r 


ROD 


R 


6 


THETA(IJ) 


THETA 


k 


RHO(IJ) 


K 


g 


(All) 




M' 


PM(1J) 




N' 


PN(1J) 




V 


PV(U) 




F 


DMAX 




k 


XLAM 


LAMBDA 


W x 


WLAM 


WEIGHT 


g\ 


DELAM 


Gl 


g2 


TCON4 


G2 


g3 


TCON5 


G3 


#4 


TCON6 


G4 


#5 


TCON7 


G5 


#6 


TCON8 


G6 


y 


Y 




kN 


XLN(I) 


LAMN 


W 


XLV(I) 


LAMV 


\ 2 M 


XLM(I) 


LAM**2*M 


kr 


RLAM 


LAMR 


kt 


TLAM 


LAMT 


(Tt 


SIGTL(I) 


SIGMA 


6 for max. oy 


JMAX 


TM 


^x/Ax 


RATIO 


WT/DEFL 


Integral in eq (11) 


YINT(IJ) 


INT 


^5 


DS 




dr 






dB 


DRDT(IJ) 




cPr 

dm 


SDRDT(IJ) 




do 


DTDS 




ds 






dr 


DRDS 




ds 
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10 


RtAD(5 
FoRMAT 


C IGO.nE.O 




RcAu(5 




READ (5 




RlAU(5 




READ(5 




KLAU(5 




RlAj(5 


20 


FORMAT 




READ (5 


30 


format 



Elastic Ring Optimization Program 

9 RUN MlTCH£*16650r02*6Q 
*Ti FOR MAlNlrMAINl 

COMMON /INTEG/F(10Q)rAREA 

DIMENSION YINT(91) 

DIMENSION PV(91) 

DIMENSION Cl(20)iC3(20)rC5(20)fC7(10)rC8(10)rXN(40)*B(91)*T(91)> 

1 k(9i)*DRDT(91)rSDRDT(9i)*RH0(91),G (91) * PM(91) f PN(91) r 

2 C2(20)rTHETA(9D 
10) Nl#N3»NSrN7rh8*NNf IGO 

(2014) 
MEANS ANOTHER SET OF DATA FOLLOWS THIS SET. 
r20) (CKI) rI=l#Nl) 
»20) (C3(I)>I=lrN3) 
r20) (C5(I) f I=lrN5) 
r20) (C7(I) f I=1*N7) 
♦20) (C8(I)*I=lrN8) 
r20) (XN(I)r I=lfNN) 
(10E8.0) 

*30) XNHUrXLLAM*S*E'GAMMA 
(5E15.0) 

c insert initial write cards here, 

C PRINT OUT ALL INPUT DATA AT START OF RUN. PROGRAM 1. 

WRITE(6#1000) NlrN3*N5'N7#N8*NN* IGO 
1000 FORMAT (19H1INTEGER PARAMETERS //(8I5)) 

WRITE (6r 1010) (CI ( I ) * 1 = 1, Nl) 

WRITE (6r 1020 )(C3( I) r 1 = 1 *N3) 

WRITE (6,1030) (C5(I)rI=lrN5) 

WRITE (6f 1040) (C7(I)rI=lfN7) 

WKITE(6,1050) (C8(I)rI=lrN8) 

WRITE(6rl060)(XN(I)rI=lrNN) 
1010 FORMAT I 3H0C1 //(8E15.5)) 
1020 FORMAT (3H0C3 //(8E15.5)) 
1030 FuRMAT (3H0C5 //(8E15.5)) 
1040 FORMAT (3H0C7 //(8E15.5)) 
1050 FORMAT (3H0C8 //(8E15.5)) 
1060 FORMAT (2H0N //(8E15.5)) 

WRITE (6#1070) 
1070 FORMATUHO HXr3H NU*l5Xr7H LAMBDA r 18X r 1HS' 19X f 1HE* 14X> 5HGAMMA ) 

WK1TE( 6rl080) XNHUrXLLAMrSrE* GAMMA 
1080 FORMAT (1H0*5E20.8) 

JCOUNT=0 

PI =3.141592654 

Do 1 I=lrNl 

C^(I)=2.0-C1(I) 



DO 


1 


J=l 


rN3 


Do 


1 


K=l 


rN5 


DO 


i 


L=l 


fN7 


DO 


1 


M=l 


rN8 


DO 


1 


N=l 


rNN 


A=Ci(I) 
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Elastic Ring Optimization Program (continued) 

D= C2<I) 
C JI GOnS FROM TO 90. 
Du 2 IJ 21*91 
Jl= iJ-1 

Tht.TA(IJ)=FLOAT(JI)*PI/180.0 
b(U)= C3(J)*(1.-C7(L)*C0S(2.*THETA(U))) 
T(Ij)= Cb(K)*(l.-C8(M)*C0S(2.*THETA(IJ))) 
TERM = A **XN(N)*SIN(THETA(IJ))**XN(N)+D **XN(N)* 

1 COS(THETA(IJ))**XN(N) 
R(U)= A*D /TERM**(1./XN(N>) 

TPi = D **XN(IM)*COS(THETA(IJ))**(XN(N)-1.0)*5IN(THETA(IJ)) 

TP2 = A **XN(N)*SIN<TNETA(IJ))**(XN(N)-1*0)*COS(THETA(IJ)) 

DKDT(IJ)= R(IJ)*(TP1-TP2)/TERM 

IF(XN(N).EQ.2.0.AND.THETA(IJ).EQ.0.0) GO TO 121 
TP3 = A ** XN(N)* 

1 <SIN(THETA(IJ))**XN<N)-(XN(N)-1.0)*SIN(THETA<IJ))** 

2 Un(N>-2.) * C0S(THETA(IJ))**2 ) 
IFCXN(N) .£Q.2.0.AND.IJ.EQ#91) GO TO 122 
GO TO 125 

121 TP3 = A**2*(-1.0) 
125 CONTINUE 

TH^= ** XN(N)* 

1 (CuS(THETA(IJ) )**XN(N)-(XN(N)-1.0)*COS(THETA(IJ) )** 

2 (XN(N)-2.>* SIN(THETA(IJ))**2 ) 
GO TO 125 

122 TP4 = D**2 * (-1.0) 
125 CONTINUE 

SDROT (IJ)= A*D * ( XN ( N ) +1 . ) * ( TP2-TP1 ) *TERM** 
1 C-2.-(l./XN(N>) )*(TP2-TPD+ A*D*TERM** (-1 •- ( 1 ./XN(N) ) ) * (TP3+TP4) 

RhO ( I J) = (R ( IJ) **2+2 . 0*DRDT ( IJ) **2-R ( IJ) *SDRDT ( IJ) ) / 
1 (R(IJ)**2+DRDT(U)**2)**(3./2.) 
C RHO IS NO* K. 

IF(T(IJ)*RH0(IJ).GE.2.0) GO TO 1 

DO 1971 JMA=lfll 

Y=-T(IJ)*.b ♦ (T(IJ)*.5) * (FLOATCJMA-D/5.0 ) 
1971 F(JMA)=Y**2/(1.0+RHO(U)*Y) 

CALL YSIMP(10.0/T(IJ)> 

YxNT(IJ)=AREA 

JXV= RHO(IJ)/T(IJ) 

DJ 1972 JMA=1»11 

Y=-T(IJ>*.5 ♦ (T(IJ)*.5) * (FL0AT(JMA-l)/5.0 ) 

SUM = (1.+ RHO(IJ)* Y) 
1972 F(JMA)=DIV**2+2.0*DIV*Y/(SUM*YINT(IJ))+Y**2/(SUM*YINT(IJ))**2 

Caul YSIMPC10.0/T(IJ)) 

G(U)= AREA 

DTub = 1./ SQRT(R(IJ)**2+DRDT(IJ)**2) 

DKUS = DRDT(1J)*DTDS 

PV(Ij)= K(IJ)*SlN(ThETA(IJ) )*DTDS-COS(THETA( IJ) )*DRDS 

Pis. ( I J) = SIN (THETA ( I J) ) *QRDS+R ( I J) *COS ( THETA ( I J) ) *DTDS 
2 CONTINUE 
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Elastic Ring Optmization Program (continued) 

C NOW WE WANT MPRIMECCALLED PM) . 
C WE DO SOME SIMPSON INTEGRATION. 

DO 101 JJ=1*91 

F(JJ)= R(JJ)*C0S(THETA(JJ))*SQRT(R(JJ)**2+DRDT(JJ)**2)*G(JJ)/B(JJ) 
1 *(1.0-XNHU**2> 

101 CONTINUE 
CALL SIMP(l.O) 
PARI =AREA 

DO 102 JJ=lf91 

F(JJ)= SIN(THETA(JJ))*DRUT(JJ)* RHO( JJ)/ (B ( JJ> *T ( JJ) ) 

102 CONTINUE 
CALL SIMP(l.O) 
PAR2 =AREA 

DO 103 JJ=lf91 

F(Jj)= COS(THETA(JJ))*R(JJ) * RHO( JJ) /(B ( JJ) *T( JJ) ) 

103 CONTINUE 
CALL SIMP(l.O) 
PAR3= AREA 

DO 104 JJ= 1*91 

F(JJ)= SQRT(R(JJ)**2 +0RUT(JJ)**2) *G(JJ)/B(JJ) 
1 *(1.Q-XNHU**2) 

104 CONTINUE 

CALL SIMP(l.O) 
PARif = AREA 
C NOW INTEGRATE TO GET W SUB LAMBDA. 
DO 105 JJ=1*91 
F(JJ)= SQRT(R(JJ)**2+DRDT(JJ)**2)*B(JJ)*T(JJ)*PI/180. 

105 CONTINUE 
CALL SIMP(l.O) 
PARb = AREA 

C USE PAR5 TO COMPUTE W SUB LAMBDA LATER ON. 
C NOW GET PM. 

DO 106 JJ=lf91 

PM(JJ)= (PAR1-PAR2-PAR3)/PAR4 -R(JJ)* COS(THETA ( JJ) ) 

106 CONTINUE 

C FORM NEW LOOP TO GET BMIN. 

JMAX=0 

DMAX =0.0 

DO 8811 IJ =1*91 

OIV = PM(IJ)*RHO(IJ)/(B(IJ)*T(IJ)) 

DENOM = DIV- PM(IJ)*T(lJ)/Cb(IJ)*(2.-RH0(lJ)*T(IJ))*YlNT(IJ)) 
1 ♦ PN(IJ)/(B(IJ)* T(IJ)) 

IF (ABS(DENOM).GT.DMAX)JMAX=IJ-l 

DivlAX = AMAX1 (ABS(UENOM)r ABS(DMAX)) 
8811 CONTINUE 

BMIN = 2.0/DMAX 

XLAM = XLLAM /(BMIN*S) 

WLAM = 4.*XLAM**2*GAMMA*PARb 
C NOW COMPUTE SOME VALUES TO COMPARE WITH CONSTANTS. 

TC0N4 = XLAM*C5(K)*(1.+ABS(C8(M))) 
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Elastic Ring Optimization Program (continued) 

ICONS = XLAM*(2.*C1(I)+T(D) 
TC0N6 = XLAM*(2.*C2(I)+T(9D) 
TCON7 = C3(J)*U.+ABS<C7<L)>) 

TCUN6 = XLAM*(2.*C2(I)-T(91) ) 
C NOW GO TO DOOPS IF CONSTRAINTS ARE NOT SATISFIED. 
JO 107 JJ=lr91 

Ti'li = (l.-XNHU**2)*G(JJ)*T(JJ)*PM(JJ)*R(JJ)*C0S(THETA(JJ) ) 
Ti«i2 = RHO(JJ)*PM(JJ)*PN(JJ) 

0S= SQRT(R(JJ)**2+DRDT(JJ)**2) 

TM3 = RHO(JJ)*PN(JJ)*R(JJ)*COS(THETA(JJ)) 

F(JJ)=1./(B(JJ)*T<JJ))*(PN(JJ)**2+12./5.*(1.+XNHU)* 
1 PV( J)**2-TM1+TM2-TM3)*DS *PI/180. 
1U7 CONTINUE 

CALL SIMP (1.0) 
AbTAR=AREA 

UtLAM=XLAM*ASTAR*BMIN*S/E 
RaT10=WLAM/0ELAM 
C COMPUTE SOME OUTPUT ITEMS. 

IF(M0D(JC0UNT»25) .NE.O) GO TO 7 
WRITE(6>4) 
«♦ F0RMAT(1H1 IX f WT/0EFL ,< +X f WEIGHT f 6X»Gl t 6X f G2 f 5X f G3 f 6X♦Gt+ f 6X t G5• 

1 7x» G6'3X» LAMBDA f 4X • A f 7X»C • 6X • D» 7X»£» 7X »F f 7X' N f 2X» TM» //) 
7 WRITE <6r 5) RATIO r WLAM'DELAM* TC0N4 r TCON5r TC0N6* TC0N7* TC0N8* XLAM* 

1 Cl(I)rC3(J) fC5(K) rC7(L) fC8(M) rXN(N) * JMAX 
5 FORMAT(lH0*2F9.2fF9.3rl2F8.3rI3) 
JCOUNT=JCOUNT+l 
1 CONTINUE 

IF(IGO.NE.O) GO TO <f4 
C WE HAVE NOW FINISHED MAJOR LOOP. 
STOP 
END 
liTl FOR SlNTrSINT 

subroutine simp(hhm) 
common /integ /f(100)#area 

H = l./(3.*HHH) 

AKEA=0.0 

OuD=o.O 

EVEN=0.0 

DO 21 I =2>90r2 

21 EVEN =EVEN ♦ F(I) 
DO 22 I =3f89»2 

22 ODD =ODO+ F(I) 

AREA = H* (F(D+F(91)+4.*EVEN+ 2.*0DD) 
RETURN 

till FOR YSUBrYSUB 

SUBROUTINE YSIMP (HHH) 
COMMON /INTEG /F(100)rAREA 
H = l./(3.*HHH) 
AREA=0.0 
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Elastic Ring Optimization Program (continued) 





OuD=0. 







EVEN=0 


.0 




DO 21 


X=2tl0t2 


21 


EVLH = 


EVEN ♦ F(I) 




DO 22 


I=3r9>2 


22 


OuO =ODD+ F(I) 




AKEA = 


H*(F(1)+F(11)+4.0*EVEN+2.0*ODD) 




RETURN 






LNO 




XUT 


MmINI 




1 


b 1 


6 3 10 


1.21 






12.7 


13. 


13.3 13.6 14. 


.305 






-.025 


0. 


.025 .05 *075 .1 


.4 


.45 


• 5 


6. 






• 30 




1000000. 150000. 30000000 


EOF 






FIN 
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Detailed Ring Analysis Program 



RUN Mi TCHEr I6b50r02rb0 
toiXT FuH MAIN2rMAIN2 

COMMON /1NTEG/FU00) f AREA 

DIMENSION XLN(9i) >XLV(9l) »XLM(9l) 

DIMENSION YINT(9i) 

DIMlNSION PVC91) 

DIMENSION 6(91) 

DIMENSION B(91) *T(9i) rR(91) * DRDT(91 ) f SDRDT (91 ) t RH0(91 ) r 
1 PM(91) f PNC91) # SIGTLC91) r THETAC9D 
10 RcAU(5r20) CI t C3* C5r C7f C8f XN 

20 Format (lOEa.o) 

IF (Ci. EG. 9999.) STOP 

RlAj(5*30) XNHUfXLLAMfSrEf GAMMA 

30 Format (seis.o) 
c insert initial print here, 
c print out all input data at start of run. program 3. 

wRlTE(6f 1000) ClrC3fC5'C7rC6rXN 
1000 FORMAT ( 1H1» 5X'2HC1 9 8Xr2HC3r 8X>2HC5* BXt 2HC1 * 8X*2HC8*8Xr lHNr // 



1 6F10.3) 
WRITE (6*1070) 

rtklTE (bflOtiO) XNHUfXLLAMrSfE^GAMMA 
: INSERT FORMATS 1070 AND 1080 HERE. 
1070 
1080 



121 
125 



.r\i rvrvi'iM i j xu i u hihu jiuou nunc . 

FuRMATUHO 10X*3H NU*l5Xr7H LAMBDA ' 18X • 1HS ' 19X f lHE'l^X* 5HGAMMA ) 

FORMAT (lH0*5E20 t 8) 

Pi =3.141592654 

C2=2 f 0-Cl 

C4 = C3*C7 

Co = C5*C8 

A =C1 

D=C2 

Do 2 I = 1*91 

JI = I-i 

THETA(I)=FLOAT(JI)*PI/l80. 

BCD = C3 * (1.-C7 *C0S(2.*THETA(I))) 

T(I) = C5 * (1.-C8 *COS(2.*THETA(I))) 

TLKM = A ** XN *SIN(THETA(I) )**XN+D **XN * 
1 COS(THETA(I))** XN 

R(I)=A*Q/(TERM **(1./XN) ) 

TP1 = D **XN* COS(THETAd) )**(XN-1.0)*SIN(THETA(D) 

TP2 = A **XN* SIN(THETA(I))**(XN-1.0)*COS(THETA(I)) 

DKDT(D=R(I) *(TP1-TP2)/TERM 

IF(XN.EQ.2.0.AND.I .EQ.l) GO TO 121 

TP3 = A **XN*(SIN(THETA(I))** XN-(XN-l.O)* SIN (THETA(I))** 
1 (XN-2.0)*C0S (THETA(I))**2 ) 

IF(XN.EQ.2*0.AND.I .EQ-91) GO TO 122 

GO TO 125 

TP3 = A**2*(-i.O) 

CONTINUE 

TP4 = D**XN *(COS(THETA(I))** XN-(XN 
1 (XN-2.0USIN (THETA(I))**2 ) 



•1.0)* COS (THETA(I))** 
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Detailed Ring Analysis Program (continued) 

Go TO 123 
1*2 TP4 = D**2 * (-1.0) 
1*3 CONTINUE 

Sl>Rl)T(I)= A*D * (XN+1.0)* (TP2-TPD* TERM** 
1 <-2.-Cl./XN))*(TP2-TPl)+ A*D*TERM**(-1.-Q./XN))*(TP3+TP<U 

RrlO(D = (R(I)**2+2.*DRDT(I)**2-R(I)*SDRDT(I)) 
1 / (R(I)**2+DRDT(I)**2)**(3. 0/2*0) 
DO 1971 J=lfll 
Y=- T ( 1 ) * • 5+ T ( I ) * • 5* ( FLOAT ( J-l ) /5 . ) 

1971 F(J)=Y**2/(1.0+RHO(I)*Y) 
CALl YSIMP(10.0/T(I)> 
YlNT(I)=AREA 

DIV = RHO(I)/T(I) 

DO 1972 J=l*ll 

Y=-T 1 1 ) * . 5+ T ( I ) * • 5* ( FLOAT ( J-l ) /5 . ) 

SUM=(1.0+RHO(1)*Y) 

1972 F(J)=0IV**2+2.*DIV*Y/(SUM*YINT(I))+Y**2/(SUM*YINT(I))**2 
CALL YSIMP(10.0/T(I)) 

G(I) = AREA 

DTDS =1./SQRT(R(I)**2 +DRDT(I)**2) 
DKDb = DRDT(I)* DTDS 

PV ( I ) = R ( I ) *SIN ( THETA ( I ) ) *DTDS-COS ( THETA ( I ) ) *DRDS 
PNli) = SIN(THETA(I))*DRDS +R ( I ) *COS(THETA (I ) ) *DTDS 
2 CONTINUE 
C NOW WE INTEGRATE TO GET MPRIME (CALLED PM) • 
DO 101 J =1*91 

F(J)=R(J)*C0S(THETA(J))*SQRT(R(J)**2+DRDT(J)**2) 
1 *G ( J) /8 ( J) * ( 1 • 0-XNHU**2 ) 

101 CONTINUE 
CALL SIMP(l.O) 
PaR1=AREA 

DO 102 J =1*91 

F(J)= SIN(THETA(J))*DRDT(J)*RHO(J)/(B(J)*T(J)) 

102 CONTINUE 
CALL SIMP(l.O) 
PAR* =AREA 

DO 103 J=lf91 

F(J)= R(J)*COS(THETA(J))*RHO(J)/(B(J)*T(J)) 

103 CONTINUE 
CALL SIMP(l.G) 
PAR3= AREA 

DO 104 J=l*91 

F(J)=SQRT(R(J)**2+DRDT(J)**2) *G( J)/B( J) * ( 1 .0-XNHU**2) 

104 CONTINUE 
CALL SIMP(l.O) 
PAR<* = AREA 

C NOW WE INTEGRATE TO GET W SUB LAMBDA. USE PAR5 LATER ON. 
Du 105 J=l*91 
F(J)=SQRT (R(J)**2 +DRDT(J)**2)* B( J> *T( J) *PI/180 . 

105 CONTINUE 
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Detailed Ring Analysis Program (continued] 

CALL SIMP(l.O) 
PAR5 = AREA 
C NOW GET PM. 

DO 106 J=l#91 

PM(J)= (PAR1-PAR2-PAR3>/PAR4 -R ( J) *COS(THETA( J) ) 

10b continue 

C FORM LOOP TO GET BMIN AND SIGTL. 
DNAX =0.0 
DO 8811 1=1*91 

DIV = PM(I)* RHO(I)/(B(I)*T(D) 

SIGTL(I)= DIV - PM(I)*T(I)/(B(I)*(2.-RH0(I)* T( I ) )*YINT( I ) ) 
1 +PN(I)/(B(I>*T(I>) 
DMAx = AMAXKABS(SlGTL(I))rABS(DMAX)) 
8611 CONTINUE 

BMIN = 2.0/OMAX 
XlAM= XLLAM /(BMIN*S ) 
FACT0R=BMIN*S*.5 
DO 6833 I =lr91 
SIGTL(I)=FACTOR*SIGTL(I) 
XLN ( I ) =XLAM*FACTOR*PN ( I ) 
XLV 1 1 ) =XLAM*FACTOR*PV ( I ) 
XLM ( I ) =XLAM**2*FACT0R*PM ( I ) 
8633 CONTINUE 

DO 107 JJ=1>91 

Thl = (l.-XNHU**2)*G(JJ)*T(JJ)*PM(JJ)*R(Jj)*C0S(THETA(JJ)> 
Tn.2 = RHO(JJ)*PM(JJ)*PN(JJ) 

DS= SQRT(R(JJ)**2+DRDT<JJ)**2) 

Th3 = RHO(JJ)*PN(JJ)*R(JJ)*COS(THETA(JJ)) 

F(JJ)=1./(B(JJ)*T(JJ))*(PN(JJ)**2+12./5.*(1.+XNHU)* 
1 PV(JJ)**2-TM1+TM2-TM3)*DS *PI/180. 
107 CONTINUE 

CALL SIMP(l.O) 
ASTAR=AREA 

DELAM=XLAM*ASTAR*BMIN*S/E 
WLAM = <*.* XLAM **2* GAMMA *PAR5 
C NOW COMPUTE SOME OUTPUT ITEMS. 
TC0|\i<4=XLAM*C5*(l.+AdS(C8) ) 
TCONb= XLAM * (2.* A +T(D) 
TC0N6= XLAM * (2.* D +T(9D) 
TCON7= C3 * (1.+ ABS(C7)) 
TCON8=XLAM*(2.0*C2-T(91) ) 
C NOW COMPUTE STRESS. FOR EACH THETA. 

C NOW WE PRINT OUT DESIRED OUTPUTr THEN READ ANOTHER DATA SET. 
WkITE(6*40) 
<*0 FORMATClHO IX f WT/DEFL f <+X • WEIGHT 1 6X'G1 • 6X »G2 f 5X f G3 f 6X'G4 f 6X f G5 f 
1 7X'G6 f 3X'LAMDA'5X*A*7X'C f 8X'D*7X , E*7X f F»7X»N f / ) 
RAT10=WLAM/DELAM 

wiKlTE<6f70)RATIOfWLAMrDELAMrTCON^fTCON5rTCON6fTCON7fTCON8fXLAMr 
1 CltC3fC5fC7rC8rXN 
70 FoRMAT(lH0r2F9.2rF9.3rl2F8.3) 



JUL. RES. SEC. C 



159 



Detailed Ring Analysis Program [continued] 



WRITE (6*80) 
60 FORMAT ( 1H0 f THETA • 6X • R f 9X f LAMR f 9X f T f 8X • LAMT f 9X • L • 9X f K • 9X • INT f 
1 7X'LAMN»7X»LAMV f 4X f LAM**2*M f 6X f SIGMA*//) 

DO 64 I =1*91 

J =1-1 

RLAM=R(I)*XLAM 

TlAm=T(I)*XLAM 

WRIT£(6*90)J*R(I) * RLAM* T( I ) * TLAM*B( I ) * RHO( I ) * YINT( I) * XLN( I ) * 
1 XLV(I)*XLM(I)*SIGTL(I) 
90 FQRMATUH I5*11E11.4) 

IF(MOD(Jr5).EQ.O.ANQ.J.NE.O) WRITE (6*110) 
64 CONTINUE 
110 FORMAT (1H ) 

GO TO 10 

STOP 

E*D 
ON FOR SINT*SINT 

SU3R0UTINE SIMP(HHH) 

COMMON /INTEG /F(100)*AREA 

H = 1.0/(3.0*HHH) 

AR£A=0.0 

OuO=0.0 

£v£N=O.Q 

00 21 I =2*90*2 
21 Ev£i\i =EVEN ♦ F(I) 

DO 22 I =3*69*2 
ZZ OjD =ODD+ F(I) 

AKtA = H* (F(1)+F(91)+4 # *EVEN+ 2.*0DD) 

RETURN. 

Ei\iO 

on for ysu6*ysub 

subroutine ysimp(hhh) 
common /integ /f(100)*area 

N = l./(3.*HHH) 

AKtA=0.0 

ODD=0.0 

EVEN=0.0 

DO 21 I =2*10*2 

21 EVEN =EVEN •► F(I) 
DO 22 I =3* 9*2 

22 ODD =ODD+ F(I) 

AREA = H* (F(D+F(11)+4.*EVEN+ 2.*0DD) 
RETURN 



£,\D 

XOT MAIN2 

1.22 13.8 

.3 

1. 14. 

• 3 

EOF 
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1000000. 
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